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ABSTRACT 

This  study  examines  the  binary  source  correlation  technique  for  deter- 
mining vertical  profiles  of  the  refractive  index  structure  parameter,  Cn^. 
Theoretical  intensity  scintillation  covariance  functions  and  power  spectra 
for  atmospheric  layers  of  depth  Ah  at  a  mean  altitude  of  h  are  derived. 
These  functions  are  related  to  the  photoelectron  counting  statistics  and  the 
spatial  covariance  of  photoelectron  counts  for  binary  point  sources.  A 
linear  detector  array  in  the  exit  pupil  of  an  optical  system  is  examined,  and 
the  effects  of  the  detection  process  uncertainty  are  explored.  A  lower  bound 
for  the  expected  error  of  an  experimental  determination  of  the  spatioangu- 
lar  covariance  of  counts  is  derived.  This  error  is  then  minimized  via  a  least 
squares  analysis  using  the  redundant  information  available  in  pairwise 
multiple  correlations  of  a  signal  from  the  detector  elements  of  a  ten  element 
linear  array.  The  refractive  index  structure  parameter  profile  is  then 
derived  and  found  to  be  the  undetermined  coefficients  of  the  spatial  covari- 
ance weighting  functions  in  the  least  squares  analysis. 
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I.  INTRODUCTION 

The  optical  remote  sensing  of  refractivity  fluctuations  induced  by 
turbulent  mixing  is  of  the  utmost  importance  to  studies  of  coherent 
electromagnetic  (EM)  radiation  propagation  in  the  atmosphere.  A  quanti- 
tative characterization  of  atmospherically  induced  refractivity  fluctuations 
as  a  function  of  altitude  is  required  in  the  design  and  system  analysis  of 
adaptive  optical  systems. 

Adaptive  optical  systems  are  systems  that  attempt  to  correct  a  perturbed 
wavefront  via  mechanical  or  electrooptic  methods.  These  methods  include 
systems  with  mechanically  deformable  optical  surfaces  and  phase 
conjugation  techniques  using  nonlinear  optical  materials.  However,  the 
implementation  of  these  techniques  to  correct  for  atmospheric  degradation 
of  optical  signal  quality  depends  on  a  thorough  probe  of  atmospheric  optical 
effects. 

Typically,  three  quantities  are  used  to  characterize  the  atmospheric 
propagation.  These  are  the  isoplanatic  angle  (60),  the  spatial  coherence 
length  of  the  atmosphere  (ro)  and  vertical  profiles  of  the  refractive  index 
structure  parameter  Cn^.  The  purpose  of  this  thesis  is  to  study  theoreti- 
cally a  single  technique,  binary  source  correlation,  with  the  intent  of  find- 
ing a  robust  optical  technique  to  profile  the  refractive  index  structure 
parameter. 

Several  different  approaches  to  the  profiling  problem  have  been 
employed  with  varying  degrees  of  success.  These  include  both  active  and 
passive  means.    Active  probes  of  atmospheric  structure  include  acoustic 


sounders  (up  to  approximately  1  km)  [Refs.  1,2]  and  pulsed  Doppler  radar 
(3-20  km)  [Refs.  3,4].  Also,  direct  profile  measurements  may  be  made  using 
microthermal  probes  mounted  on  a  balloon  [Ref.  5].  Passive  means  of  find- 
ing the  propagation  path  variations  of  the  structure  parameter  include 
direct  inversion  of  the  amplitude  scintillation  covariance  function  [Refs. 
6,7],  spatially  filtered  apertures  [Refs.  8,9]  and  the  crossed  beam  or  binary 
source  technique. 

The  crossed  beam  method  was  originally  proposed  by  Fisher  and 
Krause  [Ref.  10]  and  also  by  Wang,  Clifford  and  Ochs  [Ref.  11].  Essentially, 
the  crossed  beam  technique  proposed  that  two  optical  sources  and  two 
receivers  be  arranged  such  that  the  beams  from  the  sources  cross  at  a  point 
in  space.  The  cross  correlation  of  the  receiver  outputs  allow  the  determina- 
tion of  the  turbulence  characteristics  at  the  beam  crossing  point.  The 
crossed  beam  technique,  hereafter  known  as  the  binary  source  method,  was 
implemented  experimentally  by  several  French  teams  [Refs.  12-15]  using  a 
binary  star  for  the  source.  This  method  had  advantages  not  present  in  the 
aforementioned  passive  techniques. 

The  greatest  advantages  of  the  binary  source  technique  are  high  spatial 
resolution  and  numerical  stability  of  the  algorithm.  This  is  in  sharp 
contrast  with  techniques  involving  the  inversion  of  an  integrated  profile. 
The  inversion  of  the  scintillation  covariance  is  often  unstable  due  to  noise 
[Ref.  16];  and,  the  spatial  filtering  technique  has  limited  spatial  resolution 
due  to  the  width  of  the  altitude  weighting  functions  [Refs.  8,9].  Addition- 
ally, an  artificial  source  will  be  available  on  the  relay  mirror  experiment 
satellite  and  it  is  hoped  that  these  techniques  will  be  directly  applicable  to 
profile  determinations  made  using  time  delayed  correlations  from  this 
source.  [Ref.  17] 
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II.   THEORETICAL  BACKGROUND 

To  fully  appreciate  the  binary  source  method  requires  an  understand- 
ing of  the  theory  of  EM  propagation  in  turbulent  media.  The  seminal  work 
in  this  field  was  carried  out  in  the  Soviet  Union  in  the  1940s  and  1950s  and 
this  work  is  summarized  in  Tatarski  [Ref.  18].  Modern  treatments  of  the 
propagation  problem  in  terms  of  the  formalism  of  scattering  theory  have 
yielded  a  better  understanding  of  the  propagation  problem;  however,  this 
formalism  has  proven  impotent  in  furthering  the  predictive  power  of  the 
classical  theory  [Ref  19]. 

The  classical  theory  as  developed  by  Tatarski  and  also  summarized  by 
Clifford  [Ref  19],  relies  on  a  statistical  description  of  the  properties  of  the 
turbulent  medium.    A  summary  of  these  properties  follows. 

A.   TURBULENCE  IN  GENERAL 

Turbulence  is  a  property  of  fluid  flows  and  these  flows  are  governed  by 
the  Navier-Stokes  equations.  The  flow  is  usually  characterized  by  a  dimen- 
sionless  parameter  called  the  Reynolds  number,  Re,  with 

Re  =  vL/n.  (2.1) 

The  quantities  v,L  and  yi  are  the  characteristic  velocity,  length  scale  and 
kinematic  viscosity  respectively.  For  small  values  of  Re  viscous  dissipation 
dominates  and  the  flow  is  laminar.  As  Re  grows  in  magnitude,  then  a 
critical  value  is  reached  (approximately  1000  in  the  atmosphere)  beyond 
which  fluctuations  in  the  velocity  field  are  no  longer  damped.  For  values  of 
Re  greater  than  Re  critical  the  flow  rapidly  becomes  chaotic. 


Following  Lumely  [Ref.  20],  there  are  several  distinguishing  character- 
istics of  turbulent  flow.   These  are: 

-  irregularity 

-  large  Reynolds  numbers 

-  diffusivity 

-  dissipation 

-  three  dimensional  vorticity  fluctuations. 

The  characteristics  listed  above  are  essential  elements  of  any  turbulent 
flow.  For  example,  the  diffusivity  of  velocity  fluctuations  in  a  turbulent  flow 
and  the  rapid  mixing  inherent  in  this  process  causes  an  increase  in  the 
rates  of  momentum  and  heat  transfer. 

Turbulence  is  also  characterized  by  three  dimensional  fluctuations  of 
the  vorticity;  and,  no  flow  of  less  than  three  dimensions  is  truly  considered 
turbulent.  The  vorticity  fluctuations  take  place  at  multiple  spatial  scales 
with  each  scale  characterized  by  a  different  value  of  Re. 

The  energy  in  a  turbulent  flow  is  ultimately  dissipated  by  viscous 
friction.  This  dissipation  converts  the  macroscopic  kinematic  energy  of  the 
flow  to  thermal  energy;  and,  this  occurs  at  the  smallest  or  inner  scale  of 
turbulence.  Fluctuations  with  a  scale  dimension  less  than  the  inner  scale 
(1-10  mm  in  the  atmosphere)  are  damped  by  viscous  dissipation. 

The  largest  scale  in  a  turbulent  flow  is  determined  by  the  boundary 
conditions.  This  largest  or  outer  scale  of  turbulent  motion  is  where  energy 
is  pumped  into  the  flow.  The  energy  then  cascades  from  the  largest  to  the 
smallest  scale  sizes  adiabatically  and  the  energy  is  dissipated  at  the  inner 
scale.  The  range  of  scales  between  1©  and  Lq,  the  inner  and  outer  scales,  is 
called  the  inertial  subrange. 


The  inertial  subrange  is  so  called  because  the  time  scale  associated 
with  the  inertial  transfer  of  energy  from  larger  to  smaller  scales  is  much 
shorter  than  the  time  scale  characterizing  dissipation;  therefore,  the 
energy  transfer  is  essentially  adiabatic.  However,  the  Reynolds  number  is 
dependent  on  the  scale  and  decreases  as  the  scale  decreases.  When  a  scale 
is  reached  such  that  Re  is  less  than  the  critical  value  then  laminar  motion 
results  at  that  scale  size  and  smaller.  Using  this  fact  and  the  definition  of 
the  Reynolds  number,  the  inner  and  outer  scales  may  be  related  (within  a 
constant)  by  dimensional  analysis  [Ref.  18].   The  resulting  relationship  is 

(lo/Lo)-4/3  =  Re  (2.2) 

with  Re  the  Reynolds  number  of  the  flow  as  a  whole. 

The  existence  of  a  well  defined  inertial  subrange  provides  a  handle  by 
which  the  statistical  analysis  of  turbulence  is  carried  out.  This  approach  is 
pursued  further  in  the  following  sections. 

B.    STATISTICAL  DESCRIPTION  OF  TURBULENCE 

The  Navier-Stokes  equations  are  ill-posed  since  there  are  more  vari- 
ables than  equations;  and,  unsolvable  unless  certain  restrictive  assump- 
tions are  made.  Unfortunately,  these  same  restrictions  exclude  the 
turbulent  regime  of  solutions  and  a  statistical  approach  is  required  to 
analyze  the  flow. 

The  full  statistics  would  consist  of  a  knowledge  of  the  multidimensional 
probability  distribution  for  the  velocity  field.  However,  the  statistical  analy- 
sis is  complicated  by  the  fact  that  turbulence  is  a  nonstationary  process. 
Nonstationarity  implies  that  the  mean  and  higher  order  moments  of  the 
velocity  field  are  also  stochastic  functions  of  time.    Currently,  the  only  well 
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founded  statistical  analysis  of  the  velocity  field  statistics  is  due  to 
Kolmogorov  and  is  valid  only  for  scales  in  the  inertial  subrange  of  fully 
developed  turbulence.  [Ref  18] 

Though  the  ordinary  moments  of  the  velocity  field  distribution  are  non- 
stationary,  a  statistical  analysis  based  on  a  second  order  quantity  called  the 
structure  tensor  is  carried  out  by  Tatarski  [Ref.  18].  The  structure  tensor  is 
defined  in  general  by: 

Dij(r)  =  ([vi(r,+r)-vi(r,)]  [vj(ri+r)-vj(ri)])   ,  (2.3) 

with  the  angled  brackets  representing  an  ensemble  average.  This  quantity 
is  approximately  stationary  for  scale  sizes  in  the  inertial  subrange.  The 
structure  tensor  may  be  further  simplified  by  three  assumptions  and  these 
are: 

-  local  homogeneity 

-  local  isotropy 

-  incompressible  turbulence. 

Homogeneity  implies  that  the  structure  tensor  depends  only  on  the  dis- 
placement r  and  not  a  particular  location  in  space.  Isotropy  implies  rota- 
tional invariance  or  radial  symmetry  in  that  only  the  magnitude  of  r  and 
not  the  direction  is  important.  The  assumption  of  local  homogeneity  and 
isotropy  is  a  weaker  assumption  that  reflects  the  fact  that  the  structure 
tensor  depends  only  on  scales  very  close  to  r.  With  the  third  assumption  of 
incompressibility,  that  is  the  divergence  of  y.  equals  zero  [Ref.  21],  homo- 
geneity and  isotropy  yield  the  following  expression  for  the  structure 
function: 

Drr(r)  =  ([Vr(ri+r)-Vr(r,)]2)     ,  (2.4) 


where  Drr  is  the  projection  of  the  structure  tensor  in  the  direction  parallel 
to  the  displacement  r.  This  function  describes  the  statistics  of  homo- 
geneous, isotropic  and  incompressible  turbulence.  [Ref.  19] 

For  values  of  r  in  the  inertial  subrange  of  scales  dimensional  analysis 
of  the  structure  function  yields: 

Drr  =  Cv2r2/3  (2.5) 

with  Cv^  the  velocity  structure  parameter  [Ref.  18].    The  structure  param- 
eter is  a  measure  of  the  intensity  of  the  turbulence  and  it  is  directly  related 
to  the  refractive  index  structure  constant  mentioned  in  the  introduction. 
1.     Atmospheric  Refractivitv  Fluctuations 

To  relate  the  velocity  structure  function  to  extensive  variables  such 
as  the  temperature  or  refractivity  then  conservative  passive  additives  must 
be  considered  [Ref  18].  Conservative  passive  additives  are  those  quantities 
that  have  no  effect  on  the  statistical  analysis  of  the  turbulent  dynamics; 
however,  the  conservative  passive  additive  is  transported  and  mixed  by  the 
velocity  fluctuations.  An  example  of  a  conservative  passive  additive  is 
potential  temperature  0,  given  by: 

e  =  T(po/p)-286  (2.6) 

with  T  the  absolute  temperature  in  Kelvin,  p  the  atmospheric  pressure  and 
Po  the  pressure  at  sea  level.  This  quantity  is  independent  of  altitude. 

Following  Clifford  [Ref  19],  the  potential  temperature  fluctuation  is 
related  to  the  refractive  index  fluctuation  by: 

An  =  -79xlO-€(p/e2)Ae  (2.7) 


with  p  the  atmospheric  pressure  in  milHbars.  Since  6  is  a  conservative 
passive  additive  then  so  is  the  refractivity  n.  And  as  shown  by  Tatarski 
[Ref.  18],  conservative  passive  additive  structure  functions  obey  the  same 
power  law  as  the  velocity  statistics  except  near  the  viscous  convective 
subrange  which  is  near  the  inner  scale.  [Ref.  20]  Therefore,  the  structure 
function  for  refractive  index  fluctuations  is  given  by: 

Dn(r)  =  Cn2r2/3.  (2.8) 

The  refractive  index  structure  function  is  used  to  find  the  power  spectrum 
of  refractive  index  fluctuations  and  since  Clifford's  derivation  [Ref  19]  is 
clear  it  is  summarized  here. 

2.     Power  Spectrum  of  Refractivitv  Fluctuations 

The  refractive  index  at  a  point  r  can  be  decomposed  into  a  mean  and 
a  fluctuating  part 

n(r)  =  (n(r))  +  ni(r).  (2.9) 

If  ni  is  an  analytic  function  then  the  spatial  frequency  spectrum  can  easily 
be  determined  by  finding  the  spatial  Fourier  transform  of  ni(r).  However, 
since  ni  is  a  stochastic  function  then  the  spatial  frequency  decomposition  is 
carried  out  using  the  three  dimensional  Fourier-Stieltjes  integral 

n^(r)  =  jexp(iK»r)dN(K)  (2.10) 

Forming  the  co variance  of  the  refractive  index  fluctuations, 

B(r+r^,r^)  =  <n^(r+r^)n^(r^)>  ,  (2.11) 

the  power  spectrum  is  found  by  using  the  Wiener-Khinchin  theorem. 


The  Wiener-Khinchin  theorem  asserts  that  the  power  spectrum 
and  covariance  form  a  Fourier  transform  pair.  Using  this  and  by  invoking 
homogeneity  and  isotropy,  CHfford  [Ref.  19]  expresses  the  power  spectrum 
of  refractive  index  fluctuations  as: 


O  (K)  =  ^—  f  drrB  (r)sin(Kr)  (2.12) 

0 


This  expression  is  simplified  further  by  use  of  the  following  relation 
between  the  structure  function  and  the  covariance: 

B  (0)-B  (r)  =  ^D  (r)  (2.13) 

n  n  2      i^ 

Using  this  relationship  and  in  the  limit  that  the  outer  scale  and  inner  scale 
go  to  infinity  and  zero  respectively  the  spatial  power  spectrum  becomes: 

On(K)  =  .033Cn2K-11^3  (2.14) 

Actually,  this  expression  is  valid  only  for  spatial  wave  numbers  such  that 
27iLo~^  «  K  «  27tlo~-'^  with  the  limiting  process  used  to  derive  it  being  only 
an  analytical  convenience. 

With  these  expressions  in  hand  the  first  order  theory  of  optical 
propagation  in  the  atmosphere  is  addressed  in  the  next  section. 

C.    FIRST  ORDER  THEORY  OF  EM  PROPAGATION  IN  THE 
ATMOSPHERE 

The  first  order  theory  of  the  propagation  of  an  EM  wave  in  the  atmo- 
sphere is  essentially  equivalent  to  first  order  scattering  theory  in  which  the 
Born  approximation  allows  an  approximate  solution  to  the  wave  equation. 


This  approach  assumes  single  scattering  of  a  monochromatic  wave 
incident  on  the  atmosphere  and  it  is  valid  for  relatively  weak  turbulence. 
For  strong  turbulence  multiple  scattering  effects  cause  many  of  the  statisti- 
cal quantities  of  interest  to  saturate  to  a  constant  value;  however,  except  for 
long  propagation  paths  through  strong  turbulence  (1-2  km  in  the  late  after- 
noon) the  first  order  theory  proves  adequate  for  many  purposes. 

The  first  order  theory  assumes  that  the  atmosphere  has  unit  magnetic 
permeability  and  zero  conductivity  and  that  the  incident  EM  wave  has  a 
sinusoidal  time  dependence.  With  these  assumptions  both  Tatarski  and 
Clifford  form  the  scalar  wave  equation  for  the  propagation  of  a  component 
of  the  electric  field.  Clifford  [Ref  19]  proceeds  to  solve  this  equation  by  the 
method  of  smooth  and  small  perturbations  and  this  solution  is  summarized 
here. 

1.     Solution  of  the  Wave  Equation 

Since  the  electric  field  is  proportional  to  the  magnetic  field,  and 
assuming  that  depolarization  effects  are  negligible,  then  Clifford  shows 
that  it  is  necessary  to  analyze  only  one  component  of  the  full  vector  wave 
equation.  The  equation  to  be  solved  is: 

V2E  +  k2n2E  =  0  (2.15) 

with  E  one  of  the  components  of  the  electric  field,  k  the  wave  number  of  the 
radiation  and  n  is  the  position  dependent  atmospheric  refi*active  index. 

The  solution  to  this  equation  is  found  by  expanding  E  in  a  series  of 
decreasing  terms: 

E=Eo  +  Ei  +  E2  +  ...  (2.16) 
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with  the  zero  order  term  corresponding  to  the  unperturbed  field,  the  first 
order  term  the  single  scattering  term  and  higher  order  terms  correspond- 
ing to  multiple  scattering.  Retaining  the  terms  to  first  order  and  using 
equation  2.9  then  equation  2.15  is  resolved  into  the  two  equations: 

V2Eo  +  k2Eo  =  0  (2.17) 

V2Ei  +  k2Ei  +  2k2ni  Eq  =  0  (2.18) 

This  separation  is  accomplished  by  substituting  E  =  Eq+Ei  and  equating 
terms  of  the  same  order  to  zero  while  neglecting  terms  of  order  ni^  or 
higher.  At  this  point  it  should  be  noted  that  each  term  Em  in  the  perturba- 
tive  expansion  is  assumed  to  be  of  order  ni"^.  Following  Clifford  and 
Tatarski  [Ref.  18-19]  the  unperturbed  field  is  assumed  to  be  a  unit  ampli- 
tude plane  wave  propagating  in  the  z  direction  and  Eq  becomes: 

Eo  =  exp(ikz)    .  (2.19) 

Substituting  this  into  the  source  term  of  equation  2.18  then  this  equation 
becomes  the  nonhomogeneous  Helmholtz  equation  for  the  perturbed  electric 
field  Ei: 

V2Ei  +  k2Ei  =  -2k2nieikz  (2.20) 

The  formal  solution  of  this  equation  is  well  known  and  it  is  the  convolution 
of  the  source  term  of  equation  2.20  with  the  Helmholtz  equation  Green's 
function: 

E.(r)  =  —     dV  V—^  [2k^n^(r')e^^"']  (2.21) 

with  the  integration  over  the  scattering  volume. 
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To  simplify  the  formal  solution  further  the  Fraunhofer  approxi- 
mation is  made  in  equation  2.21.  The  Fraunhofer  approximation  assumes 
that  X  is  much  less  than  the  size  of  the  scattering  refractive  index  inhomo- 
geneity  so  that  radiation  is  scattered  in  a  small  forward  cone.  In  fact,  if  the 
varying  refractive  index  is  perceived  as  a  weak  phase  screen  then 
qualitatively  the  angle  of  scattering  will  be  a  maximum  for  inhomogeneities 
on  the  order  of  the  inner  scale.  Therefore,  from  simple  diffraction  theory 
the  maximum  scattering  angle  should  be  on  the  order  of  X/Iq  with  X,  the 
wavelength;  and,  for  inner  scales  on  the  order  of  millimeters  then  this 
angle  is  on  the  order  of  10~^  radians.  Specifically,  the  Fraunhofer  approx- 
imation assumes  that  the  vertical  distance  from  scatterer  to  receiver, 
I  z-z'  I ,  is  much  greater  than  the  transverse  displacement,  I  p-p'  I ,  from 
the  z  axis.  Making  the  Fraunhofer  approximation  means  replacing  the 
term  Ir-r'l  in  the  denominator  of  equation  2,21  by  I  z-z' I  and  expanding 
the  term  I  r-r'  I  in  a  binomial  series  retaining  terms  to  second  order  in  the 
phase.   The  resulting  expression  for  the  perturbed  field  Ei  is: 


E 
'~         2k 

vol 

The  physical  interpretation  of  this  equation  is  that  of  a  spherical  Huygen's 
wavelet  emitted  at  the  scattering  weak  phase  screen.  The  fringes  produced 
by  the  interference  of  the  Huygen's  wavelet  with  the  unscattered  plane  wave 
are  interpreted  as  the  amplitude  and  phase  perturbations  observed  in 
atmospheric  optical  propagation. 

For  the  task  at  hand  the  quantities  of  interest  are  the  fluctuations  of 
the  amplitude  and  the  intensity  of  the  optical  signal  about  the  unperturbed 
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values.  To  find  the  amplitude  fluctuations  both  Tatarski  and  Clifford  use 
the  Rytov  approximation  which  assumes  that  the  solution  of  the  wave  equa- 
tion is  of  the  form 

E  =  Aeis  ;  (2.23) 

with  A  an  amplitude  and  s  a  complex  phase.    Letting  Eq  =  Aoexp(iso)  then 

the  ratio  E/Eq  becomes: 

E 
f=l^E^  =  X-^''Pf«^oM  (2.24) 

0  0  0 

And,  taking  the  natural  log  of  equation  2.24  splits  the  ratio  into  real  and 
imaginary  parts  as  follows: 

log  ( 1  +  ^)  =  log  ( 1  +  -^)  +  i  (S-S^)  (2.25) 

0  0 

Since  Ei  is  assumed  much  less  than  Eq  and  also  that  Ai  «  Ao,  then  the 
logarithm  can  be  expanded  quite  accurately  in  a  Taylor  series.  Carrying 
out  this  expansion  and  retaining  terms  to  first  order  yields: 

E      A 

^s/  +  i(^„)  (2.26) 

0  0 

Following  Clifford  this  implies  that  the  amplitude  ratio  is  approximately 
the  real  part  of  equation  2.22  normalized  by  the  unperturbed  field  Eq.  This 
yields  the  following  expression  for  the  amplitude  ratio: 

°  vol 

Tatarski  shows  that  inherent  in  this  approximation  is  the  assumption  that 
the  amplitude  perturbation  is  small  over  the  distance  of  one  wavelength  of 
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the  radiation.  This  is  called  the  assumption  of  smooth  perturbations  and 
this  condition  is  virtually  always  satisfied  for  optical  radiation. 

Conventionally,  a  quantity  called  the  log  amplitude  is  defined  such 
that: 

x(r)  =  ln(A/Ao)  =  Ai/Aq  (2.28) 

Since  fi-om  the  weak  scattering  viewpoint  the  log  amplitude,  x(r),  is  equal  to 
the  normalized  amplitude  fluctuation  Ai/Aq  then  it  will  be  used  in  all 
further  calculations. 

To  proceed  further  the  stochastic  variable  ni(r')  in  equation  2.27  is 
expanded  in  a  two-dimensional  Fourier-Stieltjes  integral: 


n^(x')  =  j 


=  J  do(K,z')  exp  (iK»a)  (2.29) 

with  the  expansion  in  a  plane  perpendicular  to  the  propagation  direction. 
The  random  amplitude  da  now  becomes  a  function  of  z',  the  location  along 
the  propagation  path.  Substituting  this  expression  into  equation  2.27  and 
carrying  out  the  indicated  integrations,  Clifford  obtains  the  following 
expression  for  the  log  amplitude: 

x(r)  =  y^'l'(k\dzdc(E..z)  sin  [^^^^^J  )  (2.30) 

Using  this  expression  the  second  order  statistics  of  the  observed  amplitude 
fluctuations  and  the  power  spectrum  of  these  fluctuations  is  derived  in  the 
next  section. 

2.     Covariance  and  Power  Spectrvmi  of  Log  Amplitude  Fluctuations 
The  spatial  covariance  of  the  log  amplitude  fluctuations  is  found  by 
using  equation  2.30  and  the  relation: 
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B^(p,z)  =  (x(p^+p,z)  x*(pj,z))  (2.31) 

Following  Clifford  and  inserting  equation  2.30  into  equation  2.31  the  spatial 
covariance  function  is  written  as: 

B jp,z)  =  k^  jj ^^■^^^^^-<^  jj  ,,,,,,„  [ kVz)] 

2 
[^-^^]  <da(K,z')do(K,z")*)  (2.32) 


-2, 
I  iS 

sin 


with  the  superscript  *  denoting  the  complex  conjugate.  The  ensemble  aver- 
age of  the  random  complex  amplitudes  is  expressed  in  terms  of  the  two 
dimensional  spectral  density  of  the  refractive  index  fluctuations, 
Fn(K',  z'-z"),  by  using  equation  1.5  in  Tatarski  [Ref.  18]: 

(da(K,z')da(K,z")*)  =  5(K-K')F^(K',z'-z")d^Kd  V  (2.33) 

with  5(K— K')  a  three-dimensional  delta  function.  The  two-dimensional 
spectrum  of  refractivity  fluctuations  is  defined  using  the  three-dimensional 
spectrum  On  of  equation  2.14  and  is  written  as: 

oo 

F  (K',z'-z")  =       dK  <D  (K',K  )  cos[k  (z'-z")]  (2.34) 

n  —  J         z    n —      z'-z-" 

Substituting  equation  2.33  into  equation  2.32,  the  observed  spatial  covari- 
ance of  the  log  amplitude  fluctuations  becomes: 

B^(p,z)  =  J  d^K  e-'^-'^  {  k^  jl  dz'dz"  sin  [^^] 

2 

Sin  [^^^^]  F^(K,z'-z")}  (2.35) 
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This  integral  is  the  Fourier  transform  of  the  quantity  in  the  curly  brackets. 
Using  the  Wiener-Khinchin  theorem  the  expression  in  the  curly  brackets  is 
identified  with  the  power  spectrum  of  the  log  amplitude  fluctuations. 

To  proceed  further  the  experimental  geometry  as  it  bears  on  the 
evaluation  of  the  integrals  for  the  covariance  and  power  spectrum  must  be 
considered.  The  following  sections  evaluate  the  above  expressions  as  they 
apply  directly  to  the  binary  source  method  of  structure  constant  profiling. 


16 


III.  THE  ANALYSIS  OF  THE  EXPERIMENT 

The  purpose  of  this  calculation  is  to  develop  the  binary  source  technique 
theoretically  for  application  to  a  linear  array  of  photosensitive  or  photo- 
emissive  detectors  in  the  exit  pupil  of  an  optical  imaging  system.  The 
optical  system  is  assumed  to  be  well  corrected  with  an  approximately  con- 
stant modulation  transfer  function  over  the  range  of  spatial  wave  numbers 
included  in  the  inertial  subrange  of  turbulence.  With  this  assumption,  the 
intensity  distribution  in  the  exit  pupil  is  related  to  the  intensity  distribution 
in  the  aperture  by  a  simple  scaling.  Therefore,  this  calculation  examines 
the  experimental  geometry  illustrated  in  Figure  1.  The  exact  scaling  rela- 
tions and  corrections  for  a  nonideal  modulation  transfer  function  will  be 
developed  in  a  later  section. 

A.   GEOMETRICAL  CONSIDERATIONS 

The  binary  source  technique  determines  the  Cn^  profile  by  finding  the 
covariance  of  scintillations  due  to  a  finite  number  of  atmospheric  layers  Ah. 
And,  a  consideration  of  the  geometry  of  Figure  1  is  used  to  simplify  the 
analysis. 

The  prototypical  binary  source  of  Figure  1  is  a  binary  star  with  compon- 
ents labeled  1  and  2.  Two  portions  of  the  wave  fronts  intercepted  by  the 
detector  elements  are  illustrated  and  the  following  quantities  are  defined: 

d  =  11-12  with  d.  the  detector  element  separation  from  center  to  center. 

Ad,  the  detector  element  width. 

h,  the  mean  crossing  altitude  of  the  wavefronts  as  illustrated. 

Ah,  width  of  crossing. 
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h2,hi;  the  maximum  and  minimum  altitudes  of  crossing,  respectively. 
0,  the  angular  separation  of  the  binary  components  1  and  2. 


2       1      BINARY 
SOURCE 


ALTITUDE 


Figure  1.   The  Experimental  Geometry 
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With  these  definitions  and  examining  the  geometry  of  Figure  1,  two 
relationships  between  these  variables  are  apparent: 

d  =  eh  (3.1) 

Ad  =  eAh/2  (3.2) 

Equation  3.1  gives  the  connection  between  the  detector  element  separation, 
d,  and  the  crossing  altitude,  h.  The  second  equation  gives  the  relation 
between  the  detector  element  width  and  the  width  of  the  crossing  at  a  mean 
height  of  h.  Equation  3.2  is  simply  derived  by  considering  the  difference 
between  the  maximum  and  minimum  crossing  altitudes  and  using  equa- 
tion 3.1.  Note  that  equation  3.1  assumes  that  the  detector  array  is  parallel  to 
a  line  segment  joining  the  apparent  positions  of  the  binary  source  compo- 
nents.  If  this  condition  is  not  satisfied  then  equation  3.1  becomes: 

d  =  eh  sec  (|)  (3.3) 

with  <t)  the  azimuthal  angle  between  the  array  and  the  line  joining  the 
source  components.  Equation  3.1  is  assumed  valid  for  this  calculation; 
since,  the  optical  system  can  in  principle  be  aligned  and  driven  such  that 
this  relation  holds. 

B.    POWER  SPECTRUM  OF  LOG  AMPLITUDE  FLUCTUATIONS  IN  Ah 

Given  a  picture  of  the  experimental  geometry  the  derivation  of  the 
power  spectrum  of  log  amplitude  fluctuations  for  a  layer  Ah  at  an  altitude  of 

h  is  carried  out  using  equation  2.35.  This  power  spectrum  is  then  used  to 
find  the  expected  spatial  covariance  of  the  log  amplitude  fluctuations  and 
the  related  covariance  of  intensity  scintillations  via  the  Wiener-Khinchin 
theorem. 
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For  a  layer  Ah,  isotropy  and  equation  2,35  imply  the  following  expres- 
sion for  the  power  spectrum  of  the  log  amplitude  fluctuations: 


Ah 


F^(K,0)  =  k^JJ. 


dz'dz"  sin  K 


2(z-z') 
2k 


•     K  (z-z")  , 

sin  — F^(K,z  -z  ) 

2k  n 


(3.4) 


with  Fn(K,z'-z")  equal  to  Fn(K,z"-z')  the  two-dimensional  power  spectrum 
of  refractive  index  fluctuations  [Ref.  19].  Also  note  that  the  origin  of  the 
coordinate  system  is  at  h2  (z=0)  and  the  integration  is  over  the  scattering 
volume  or  layer. 

To  make  full  use  of  the  symmetry  in  equation  3.4  the  following  variable 
change  is  made  [Ref.  18]: 


^  =  z'-z"    and    2ri  =  z'-i-z" 


(3.5) 


Using  the  Jacobian  of  this  transformation  then, 


dz'dz"  = 


dzX   dz/ 


y2 

-1/2 

d^dn  = 

1 

1 

d^dn 


(3.6) 


therefore,  dz'dz"  =  d^dr).    Using  this  variable  change  and  the  trigonometric 
identity  sin(a-J-b)sin(a-b)  =  sin^a-sin^b  with, 


a  =  ^z-Ti);      b  =  —  m 


2k 


(3.7) 


then  equation  3.4  becomes: 
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^^(K,0)  =  k^J  J  d^diiF^(K,^)[sin2(-^-^  )-sin^(^)]         ^^.S) 


The  integrand  is  simplified  further  by  noting  Tatarski's  observation  on 
page  135  of  Reference  18  that  Fn(K,^)  rapidly  approaches  zero  for  K^  greater 
than  or  equal  to  unity.  This  derives  fi*om  the  fact  that  fluctuations  in  the 
refractive  index  separated  by  displacements  greater  than  Lq,  the  outer 
scale,  are  uncorrelated;  and,  by  the  Wiener-Khinchin  theorem  the  power 
spectrimi  vanishes  as  the  corresponding  covariance  of  fluctuations  goes  to 
zero.  Also  note  that  the  assumption  of  smooth  perturbations,  X  «  Iq, 
implies  that  k  »  K  over  the  significant  part  of  the  integrand  of  equation  3.8. 
Therefore,  over  the  significant  part  of  the  range  of  integration  then: 

4r^«l         and        sin^^-0        .  (3.9) 

4k  4k 

And,  equation  3.8  simplifies  to: 

F^(K,0)  =  k^  J  J  d^dri  F^(K,^)  sin^  (^^-^^)  (3-10) 

The  evaluation  of  this  integral  depends  on  the  assumption  of  a  form  for 
Fn(K,^)  and  this  form  depends  on  the  magnitude  of  Ah  compared  with  the 
inner  and  outer  scale.  For  the  experimental  geometry  of  Figure  1  the  ten 
detector  elements  of  the  linear  array  will  divide  the  total  atmospheric 
volume  into  ten  separate  non-overlapping  regions.  The  width  of  these 
regions  will  be  on  the  order  of  kilometers  for  observations  of  the  profile  from 
the  planetary  boundary  layer  to  about  20  kilometers.  This  is  much  greater 
than  outer  scales  observed  in  the  atmosphere  [Ref.  2,19]  and  it  is  safely 
assumed  that  Ah  »  Lq. 
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Under  the  conditions  discussed  in  the  previous  paragraph  the  layer  Ah 
includes  many  separate,  uncorrelated  regions  of  turbulence.  Following 
Tatarski  [Ref.  18]  the  power  spectrum  of  refractive  index  fluctuations  for 
sublayers  of  Ah  on  the  order  of  the  outer  scale  is  still  given  by  equation  2.14; 
however,  the  refractive  index  structure  constant  is  now  a  function  of  r\. 
Therefore,  the  refractive  index  fluctuation  power  spectrum  takes  the  form: 


n    '     n 


F^(K,tTi)  =  C^(7i)  f JK,^)  (3.11) 

and,  it  should  be  noted  that: 


J 


f  (K,^)d^  =  7i(.033)  K  ^^^=  7C  <D°  (K)  (3.12) 

n  "^  n 


with  this  expression  completely  analogous  to  equation  2.14.  With  an 
expression  for  the  refractive  index  power  spectrum  available  further 
progress  in  evaluating  the  integral  of  equation  3.10  is  possible. 

Noting  that  equation  3.10  is  even  in  ^  and  symmetric  in  r[  then  substi- 
tuting equation  3.11  into  equation  3.10  and  employing  symmetry  yields: 

Ah/2  2n 

F^(K,0)  =  2k'  J  dTiC^(Ti)  sin'  (^^-^)  J  f,(K,^)d^ 

O  0 

Ah  2(Ah-Ti) 

+  2k'  I  dr|(f  (Ti)  sin'  {^^^)  J  f,(K,^)d^  (3.13) 

Ah/2  o 

for  the  log  amplitude  power  spectrimi.  Over  most  of  the  region  of  integra- 
tion both  2r|  and  2(Ah-r|)  have  magnitudes  between  the  outer  scale  and  Ah. 
As    previously   noted,    the    refractive   index    power    spectrum    rapidly 
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approaches  zero  for  K^  greater  than  or  equal  to  one  [Ref.  18],  therefore  the 
upper  Hmits  for  the  ^  integration  are  replaced  with  infinity  to  a  high  degree 
of  accuracy.  Replacing  the  upper  limits  of  the  £,  integrals  of  equation  3.13 
and  carrying  out  the  integrations  using  equation  3.12  yields  the  following 
expression  for  the  power  spectrum: 

Ah 

F^(K,0)  =  2k\j 0° (K)C^(Ti)  sin^[^^"^^]dTi  (3.14) 

0 

Since  the  detector  elements  have  a  finite  width  then  by  equation  3.2 
spatial  variations  of  Cn^  smaller  than  Ah  cannot  be  resolved.  In  effect  the 
detector  observes  an  average  structure  constant  for  the  layer.  If  the  Cn^ 
profile  is  known  apriori  then  this  average  can  be  carried  out  in  the  Fourier 
transform  or  spatial  covariance  domain  as  an  aperture  average  over  the 
finite  detector  element  area.  However,  since  a  measurement  of  the  struc- 
ture parameter  profile  is  the  objective  then  it  is  the  power  spectrum  that 
must  be  averaged  over  the  layer  depth.  This  average  of  equation  3,14  is 
indicated  as: 

Ah 
(f^(K,0)>  =  2k\  Jo;(K)  {c^^(Ti)  sin^L^^P"^^)^"^  ^^'^^^ 


with  the  angled  brackets  denoting  the  spatial  average  over  the  layer.  The  r; 
integration  is  trivially: 

o 
(f^(K,0)>  =  2k\  0°(K)Ah  (c^(Ti)  sin^[^^^^]>  (3.16) 
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Past  observations  of  Cn^  indicate  that  it  fluctuates  rapidly  in  space  with  a 
mean  that  decreases  slowly  with  altitude  [Refs.  2,19].  Since  the  spatial 
fluctuations  of  Cn^  appear  uncorrelated  with  the  sin^  weighting  function 
then  the  average  indicated  in  equation  3.16  is  rewritten  as: 

(f^(K,0)>  =  2k\  o;(K)Ah  (c^(Ti)>{sin2 [^^^^] )  ^3.17) 

Carrying  out  the  spatial  average  explicitly  for  the  sin^  weighting  function 

yields: 

Ah 

(f^(K,0)>  =  2k\AhO°  (K)(c^(Ti))—  I  sin^  [^^^]  d^^         (3.18) 

0 

or,  for  z,  the  observer  coordinate,  equals  h2: 

{f^(K,0))  =  2k^7u  AhcD°  (K)  {c^(Tl)) 


-cosi^- 
AhK^ 


•[i-^-(¥)=^"(^)] 


with  (Cn^CTl))  the  average  or  effective  value  of  the  structure  parameter  for 
the  layer  at  a  mean  altitude  of  h.  This  is  the  observed  average  power 
spectrum  of  the  log  amplitude  fluctuations  due  to  turbulent  mixing  in  a 
layer  Ah. 

The  power  spectrum  used  by  the  various  French  teams  [Refs.  12-15]  is 
in  this  calculation's  notation  written  as: 

F  (K,0)  =  2k^7uAhO°(K)(c^(Ti))sin^(^)  (3.20) 

X  n         ^    n       '  \  2k  ^ 
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Approximating  the  sin(K2Ah/2k)  term  in  equation  3.19  by  AhK2/2k  (k  »  K) 
then  3.19  becomes: 

2k 


(f^(K,0))  =  2k\AhO°(K)(c^(Ti))  sin^(^)  (3.21) 


an  expression  identical  to  the  French  spectrum.  The  two  spectra  are  illus- 
trated in  Figure  2  and  it  is  evident  that  the  French  spectrum  is  an  excellent 
approximation  to  the  spectrum  of  equation  3.19.  Therefore,  the  spectrum  of 
equation  3.21  will  be  used  for  further  calculations  since  it  is  analytically 
simpler. 

C.    THE  SPATIOANGULAR  CORRELATION  OF  INTENSITY 
FLUCTUATIONS 

The  power  spectral  density  of  equation  3.21  can  be  Fourier  transformed 
to  yield  the  spatial  covariance  of  the  log  amplitude  fluctuations  of  a  single 
point  source  due  to  a  layer  Ah.  The  solution  of  the  wave  equation  in  Section 
II  assumed  a  plane  wave  solution  and  as  is  well  known  from  classical 
optics  a  point  source  at  infinity  produces  plane  waves  at  the  detector.  For 
the  stellar  sources  of  Figure  1,  the  transverse  spatial  coherence  function  is 
essentially  unity  for  displacements  of  less  than  five  meters  as  shown  by 
Hanbury-Brown  and  Twiss  for  Sirius  [Ref.  22].  Therefore,  an  assumption 
of  initial  spatial  coherence  of  the  source  is  excellent  for  spatial  separations 
of  less  than  five  meters.  However,  the  log  amplitude  covariance  is  not 
directly  obtained  by  cross-correlation  of  the  signal  from  the  detector 
elements. 

For  a  real  detector  operating  at  optical  fi"equencies  the  detector  response 
is  typically  proportional  to  the  intensity  and  it  is  the  covariance  of  intensity 
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Figure  2.     The  French  Spectrum  (x),  and  the  Spectrum  of 
Equation  3.19  (-) 


26 


fluctuations  that  is  probed  by  the  detector  system.  The  details  of  the 
detection  process  will  be  treated  in  later  sections;  however,  the  intensity 
scintillations  in  the  aperture  are  independent  of  the  detection  process.  The 
spatioangular  covariance  of  intensity  fluctuations  due  to  the  binary  stellar 
source  will  now  be  derived  and  related  to  the  log  amplitude  statistics. 
1.     The  Spatial  Covariance  Functions 

For  calculational  simplicity,  consider  the  wavefronts  and  two 
detector  elements  of  Figure  3.  The  normalized  spatial  covariance  of  the 
intensity  is  given  by: 

<[l(rp-<I(r  ))][l(r^)-<I(r^))]) 
Cj(p)  = (3.22) 

<I(r,))  (id^)) 

with  the  angled  brackets  denoting  an  ensemble  average  and  Kii)  the 
instantaneous  intensity  at  n-  The  sources  labeled  1  and  2  are  assumed 
independent;  therefore, 

Klj)  =  Iidj)  +  l2(rj)  (3.23) 

with  the  subscripts  1  and  2  denoting  the  sources.  Note  that  conservation  of 
energy  requires  that  the  long  term  average  illumination  in  the  exit  pupil  is 
uniform.  Substituting  this  expression  into  equation  3.22  and  rearranging 
terms  yields: 


+  (si^d^  5l2(r2^)  ^  (sii(l2)SVli))]  (3.24) 
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Figure  3.    The  Wavefronts  Intercepted  by  Detector  Elements 
at  ri,  12  Due  to  the  Binary  Source 
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with  5I=I-(I>.   Defining  the  relative  irradiance  fluctuations  as: 

5I.(r,) 
l.(r.)  =       '^  (3.25) 

and  substituting  into  equation  3.24  gives  the  following  expression  for  the 
normalized  spatial  covariance  of  intensity  scintillations: 

((Ij)  +  <l2»  ((ll)  +  (12)) 

^  ..    X  '  !    V.2   Kli(ri)  l^Cl^))  +  (li(l2)  l2<Ei))]  (3-26) 

((ij) + (I2)) 

The  first  two  terms  of  equation  3.26  are  the  weighted  autocovariances  of 
so'urces  1  and  2  with  the  weighting  functions  carr5dng  the  brightness  or 
magnitude  difference  between  the  two  stellar  components.  Since  assuming 
homogeneity  implies  that  the  normalized  autocovariances  are  equal  then 
equation  3.26  is  further  simplified  to: 

[{\(L^)  l^ir^))  +  {\(r^)  yii))]  (3.27) 

The  first  term  of  this  equation  is  the  normalized  autocovariance  of  a  single 
source  and  the  second  term  is  the  cross  covariances  of  source  components 
1  and  2. 
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The  intensity  covariances  in  the  observation  plane  are  related  to  the 
log  amplitude  covariances  by  a  relation  due  to  Fried  [Ref.  23],  and  since  his 
derivation  is  clear  and  readily  available  it  will  not  be  repeated  here.  In  this 
calculation's  notation  Fried's  result  is: 

Ci(p)  =  exp[4Bx(p)]  -  1  (3.28) 

Since  the  perturbative  solution  of  the  wave  equation  assumes  weak  turbu- 
lence and  this  implies  that  the  log  amplitude  fluctuations  are  small  then 
equation  3.28  is  expanded  in  a  Maclaurin  series  to  first  order.   The  result  is: 

Ci(p)  =  4Bx(p)  (3.29) 

Both  expressions  for  the  spatial  covariance  of  intensity  are  graphed 
in  Figure  4  as  a  function  of  the  log  amplitude  covariance.  The  first  order 
approximation  is  selected  for  further  calculations  based  on  comparison 
with  the  experimental  curve  of  Figure  4.  As  is  evident  from  Figure  4 
experimental  evidence  indicates  that  the  intensity  covariance  saturates  to 
unity  as  the  intensity  of  turbulence  increases.  Therefore,  the  first  order 
approximation  of  equation  3.29  is  actually  a  better  approximation  over  the 
range  of  validity  of  the  single  scattering  theory.  [Ref.  24] 

Using  equation  3.29  and  examining  the  terms  of  equation  3.27  indi- 
vidually then  the  spatial  covariance  of  intensity  scintillations  can  be  cast  in 
a  more  transparent  form  in  terms  of  the  log  amplitude  covariance.  For  the 
first  or  autocovariance  term  of  equation  3.27,  direct  substitution  of  equation 
3.29  yields: 

(l^d^)  \(t^))  =  4B^^(p)  =  4B^(p)  (3.30) 

However,  the  cross  covariance  terms  are  not  treated  as  straightforwardly. 
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Figure  4.     The  Normalized  Variance  of  Intensity  Scintillations  Versus 
the  Log  Amplitude  Variance  for  Various  Forms  for  Cj,  and 

Experimental  Data  Points  (x) 
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The  first  cross  term  in  equation  3.27  is  expressed  as  a  function  of 
the  log  amphtude  using  equation  3.29: 

{\(t^\(t^))  =  4B^      (p)  (3.31) 

2, 3 

This  term  is  the  cross  covariance  of  the  light  propagating  along  the  paths 
labeled  by  the  numbers  2  and  3  in  Figure  3.  Using  equations  2.31  to  2.33  and 
the  assumption  of  statistical  homogeneity  then  the  cross  covariance  of  log 
amplitude  fluctuations  is  expressed  as: 

B^     (p)  =  J  exp  [iK-rd^-r^  +  p]]  F^(K,0)d^K  (3.32) 

2,  3 

-oo 

The  assumption  of  statistical  homogeneity  in  the  layer  Ah  implies  that 
fluctuations  in  xi  are  mirrored  by  fluctuations  in  X2  at  the  crossing. 
Further  simplification  of  equation  3.32  is  achieved  by  using  equation  3.1. 
Substituting  this  relation  between  the  mean  altitude  and  the  detector 
separation  in  equation  3.32,  this  cross  term  becomes  the  same  function  as 
the  autocorrelation  but  with  the  origin  shifted  by  9h: 

B       (p)  =  B  (p-0h) ;      r_-r,  =  -  eh  (3.33) 

X23  X  z      1 

An  exactly  analogous  derivation  for  the  second  cross  term  yields: 

B^     (p)  =  B^(p+eh) ;     1^-12  =  eh  (3.34) 

1.4 

That  is  light  propagating  along  paths  1  £ind  4  in  Figure  3  have  a  maximum 
in  the  spatial  covariance  for  refractive  index  inhomogeneities  with  a  scale 

of  eh. 
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Substituting  equations  3.33  and  3.34  in  equation  3.27  then  the  spatial 
covariance  of  intensity  scintillations  is  expressed  as: 

Cl(P>  =       ^  4B  (p)  + i ^ 

•  4[B^(p-eh)  +  B^(p+eh)]  (3.35) 

with  the  log  amplitude  covariance  given  by  the  two  dimensional  Fourier 
transform  of  the  power  spectrum  of  equation  3.21.    Also  note  that  equation 
3.35  depends  explicitly  on  0  and  the  covariance  is  now  called  the  spatio- 
angular  covariance  as  a  remainder  of  this  dependence. 
2.     Evaluation  of  the  Spatioangrular  Covariances 

The  normalized  spatial  covariance  of  equation  3.35  indicates  that 
this  intensity  expression  can  be  constructed  from  the  log  amplitude  auto- 
covariance  by  appropriate  scalings  and  corresponding  shifts  of  the  origin. 

Using  the  Wiener-Khinchin  theorem  the  log  amplitude  autocovari- 
ance  in  the  aperture  plane  is  given  by: 

oo 

B  (p)  =    expCiK^a]  F  (K,0)d^K  (3.36) 

A.  J  A 

— oo 

with  the  power  spectrum  given  by  equation  3.21,  and  letting  the  average 
spectrum  condition  be  understood.  Note  that  {CnHr[))  is  now  identified  as 
the  mean  or  effective  structure  parameter  in  the  layer  at  a  mean  altitude  of 
h.  Equation  3.36  is  the  two-dimensional  Fourier  transform  of  the  power 
spectnmi.  However,  the  assumption  of  isotropy  and  an  assumption  of  a 
circular  aperture  imply  that  equation  3.36  can  be  recast  as  a  Hankel 
transform: 
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circular  aperture  imply  that  equation  3.36  can  be  recast  as  a  Hankel 
transform: 


oo 

=  27U  Jf^(K, 


B  (p)  =  27U  F  (K,0)J  (Kp)KdK  (3.37) 


0 
0 


with  Jo(  )  the  zeroth  order  Bessel  function.  Substituting  from  equation  3.21, 
the  autocovariance  becomes: 

oo 

B^(p)  =  Mjj^(Kp)K"^^sin^(^)dK  (3.38) 

0 

with  M=47i:2(.033)k2Ah<Cn2>.  Note  that  this  integral  is  vahd  only  for  the 
inertial  subrange  of  spatial  scales  since  the  power  spectrum  is  divergent  at 
the  origin.  Therefore,  the  limits  of  integration  are  changed  to  K'=27c/Lo  and 
k"=27C/lo,  with  Iq  and  Lq  the  inner  and  outer  scales  respectively. 

In  the  limit  that  p  goes  to  zero  then  equation  3.38  becomes: 

K" 

B^(o)  =  M  J  K'^^^sin^  i~w)  ^  (3.39) 

K' 

and  this  expression  is  the  variance  of  the  log  amplitude  fluctuations.    If 


yXh,  with  X  the  EM  wavelength,  is  much  greater  than  unity  then  the  limits 
of  integration  may  be  extended  from  zero  to  infinity  with  little  error 
[Ref  18].  Equation  3.39  becomes: 

oo 

B^(0)  =  mJ K"^^^sin^(^)  dK  (3.40) 
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and  this  expression  has  been  evaluated  by  Tatarski  in  Chapter  8  of  Refer- 
ence 18.  In  the  notation  of  this  calculation  the  variance  of  log  amplitude 
fluctuations  for  the  layer  at  a  mean  altitude  of  h  is  given  by: 

B  (0)  =  .564  Ah  h^''^  k*^^^  ((?  )  (3.41) 

The  integral  of  equation  3.38  is  not  analytic  and  a  numerical  inte- 
gration must  be  carried  out  to  find  the  spatial  covariance.  Normalizing 
equation  3.38  by  the  variance  and  carrjring  out  a  numerical  integration 
yields  the  function  of  Figure  5  with  p  measured  in  units  of  V?ih.  [Ref.  19] 

The  intensity  spatial  covariance  is  constructed  by  appropriate  shift- 
ing and  scaling  of  the  log  amplitude  covariance  as  indicated  by  equation 
3.35.  Assuming  that  the  components  of  the  binary  source  are  of  equal 
intensity  then  equation  3.35  becomes: 

Cj(p)  =  2B^(p)  +  [B^(p-0h)  +  B^(p+0h)]  (3.42) 

Normalizing  each  term  of  this  equation  separately: 

Cj(p)  =  2b^(p)  +  b^(p+0h)  +  b^(p-0h)  (3.43) 

with  bx(p±0h)=Bx(p±0h)/Bx(O).  A  graph  of  this  function  for  h  of  ten  kilo- 
meters is  constructed  in  Figure  6  with  p,h  scaled  by  VXh. 

For  experimental  implementation  the  needed  covariance  functions 
could  be  computed  beforehand  and  stored  for  later  use.  Another  alternative 
would  be  to  find  a  suitable  Chebyshev  polynomial  expansion  of  the  integral 
of  equation  3.38.  However,  the  effects  of  nonideal  optics  must  be  considered 
and  this  implies  a  numerical  calculation  of  the  log  amplitude 
autocovariance. 
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Figure  5.    The  Normalized  Log  Amplitude  Spatial  Covariance 
of  Fluctuations  [Ref.  19] 
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Figure  6.    The  Normalized  Spatial  Covariance  of  Intensity 
Scintillations 
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3.     Nonideal  Optical  Corrections  to  the  Autocovariance 

The  actual  detector  system  is  coupled  to  an  optical  system  with  the 
real  detector  array  in  the  exit  pupil.  Taking  the  linear  systems  viewpoint  of 
the  optical  system  then  it  is  viewed  as  a  linear  transformation  of  the  optical 
input  at  the  aperture.  The  individual  spatial  frequencies  in  the  aperture 
are  transformed  to  the  exit  pupil  via  the  optical  transfer  function  (OTF). 
The  OTF  of  the  optical  system  is  defined  as  the  ratio  of  the  output  spatial 
power  spectrum  to  the  input  spectrum,  and  it  is  a  spatial  frequency  depen- 
dent complex  quantity  with  the  modulus  defining  the  modulation  transfer 
function  (MTF)  and  phase,  the  phase  transfer  function  (PTF).  The  optical 
system  is  assumed  on  axis  and  since  phase  shifts  in  centered  optical 
systems  occur  off  axis  then  the  MTF  is  the  quantity  of  interest.  The  MTF 
describes  the  filtering  of  spatial  frequencies  by  the  optical  system.  For  the 
detector  array  in  the  exit  pupil  a  spatially  modulated  input  will  result  in  a 
spatially  modulated  output  that  is  suitably  scaled;  however,  if  the  MTF  falls 
to  zero  (optical  system  cutoff)  rapidly  or  is  a  rapidly  var5dng  function  of 
spatial  fi-equency  then  the  output  spatial  spectrum  is  found  from 


F  (K,0)  I     ,    ,  =  MTF(K)  F  (K,0)  I  .      ,  (3.44) 

X  output  X  input 


by  the  definition  of  the  OTF.    The  spatial  autocovariance  then  becomes  the 
Hankel  transform  of  equation  3.44  as  indicated  in  the  previous  section. 

The  power  spectrum  of  log  amplitude  fluctuations  in  the  aperture  is 
illustrated  in  Figure  2  and  it  is  noted  that  this  spectrum  falls  rapidly  to  zero 
as  a  function  of  spatial  fi-equency.  Therefore,  the  assumption  of  unity  MTF 
is  quite  good  for  reasonable  aperture  sizes  (greater  than  approximately 
20  cm).     The  MTF  will  be  assumed  constant  for  the  balance  of  this 
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calculation;  however,  experimental  implementation  of  this  scheme  should 
include  a  measurement  of  the  MTF  of  the  optical  system. 

With  the  assumption  of  constant  MTF  over  the  inertial  subrange, 
the  scaling  between  the  aperture  and  exit  pupil  is  derived  from  geometrical 
considerations.  As  is  apparent  from  Figure  7  the  average  intensity  in  the 
exit  pupil  is  increased  by  the  factor: 

tX  I.  (3.45) 


'1^. 


E 

with  Da  the  aperture  diameter,  De  the  diameter  of  the  exit  pupil  and  Tr  the 
overall  transmittance  of  the  optical  system.  For  the  linear  detector  array 
illustrated,  the  scaled  array  in  the  aperture  is  characterized  by: 

A.=A,d2,  I  ,      ;    p,  =  p^D,/  (3.46) 


with  A  the  detector  element  area,  p  a  linear  displacement  along  the  array 
and  the  subscripts  A  and  E  referring  to  the  aperture  and  exit  pupil  respec- 
tively. For  the  remainder  of  this  calculation  the  optical  system  coupled  to 
the  detector  is  assumed  to  have  unity  transmittance  and  constant  MTF  over 
the  inertial  subrange.  Under  these  assumptions  the  covariances  of  the 
previous  section  are  assumed  valid. 

Given  the  assumption  of  an  ideal  optical  system  the  interaction  of 
the  electromagnetic  field  with  the  detector  elements  must  now  be  consid- 
ered. The  next  section  discusses  the  photon  counting  statistics  of  the  binary 
source  profiling  method. 
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IV.  PHOTOELECTRON  COUNTING  STATISTICS 

The  intensity  scintillations  are  detected  by  a  linear  array  of  photosensi- 
tive or  photoemissive  detectors,  and  the  discrete  nature  of  this  detection 
process  must  be  considered.  The  incident  electromagnetic  flux  causes  the 
photosensitive  surface  of  a  detector  element  to  eject  a  photoelectron  and 
these  electrons  are  then  collected  and  analyzed  as  the  detector  element 
signal.  Following  Saleh  [Ref.  25],  this  derivation  treats  the  electromagnetic 
field  classically  as  a  stochastic  quantity;  and,  the  photoelectron  production 
process  is  treated  semiclassically.  This  approach  does  not  invoke  the 
photon  concept  and  has  the  advantage  of  not  demanding  the  full  quantum 
electrodynamic  treatment  of  the  electromagnetic  field.  The  semiclassical 
approach  is  compatible  with  the  full  quantum  analysis  of  the  photoemission 
process  provided  certain  mildly  restrictive  assimiptions  hold  true.  [Ref.  25] 

To  proceed  with  this  analysis  a  certain  amount  of  statistical  back- 
ground is  required.  The  definitions  and  relations  used  in  this  calculation 
are  summarized  in  the  next  section  and  the  notation  used  is  that  of 
Reference  25. 

A.   STATISTICAL  BACKGROUND 

For  an  experiment  that  counts  photoelectrons  the  observed  data  is  the 
number  of  counts  (ni,n2,...nni)  in  the  intervals  (tj,tj-i-T)  with  J  equal  to 
(1,2,. ..m).  Since  the  photoemission  process  is  assumed  probabilistic,  the 
number  of  counts  in  disjoint  time  intervals  are  assumed  independent. 
Therefore,  a  rate  density  P(t)  is  defined: 
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(n  (t  +  At)-N.(t)) 

p(t)  =  lim   — (4.1) 

At^O  At 

with  NxCt)  the  number  of  events  in  the  interval  T.    This  function  is  also 
known  as  the  arrival  rate  of  a  Poisson  process  [Ref.  26]. 

Given  an  assumption  of  independence  in  disjoint  intervals  then  the 
process  described  is  called  a  Poisson  point  process  (P.PP)  with  rate  density 
P(t).  The  number  of  points  (counts)  in  a  single  interval  (t,t+T)  has  a  Poisson 
density: 

P(n)  =  -V  exp(-W)  (4.2) 

n! 

with  W,  the  integrated  rate  of  events  defined  by  the  relation: 

t  +T 

0 


J 


W  =  J  p(t)dt  (4.3) 

f 

0 

And,  the  time  T  is  the  counting  time.    If  (3  is  a  constant  then  the  P.PP  is 
called  homogeneous  and  equation  4.2  becomes: 

P(n)  =  -^£^exp(-pT)  (4.4) 


n 


The  integer  random  variable  n  is  also  described  statistically  by  the 
moment  generating  function  (mgf).  The  ordinary  moment  generating 
function  is  defined  by: 

oo 

Q  (s)  =  (exp  (-ns))  =  Z       exp  (-ns)  P(n)  (4.5) 

n  ^  '  n=o 

with  P(n)  given  by  equation  4.2.    Carrying  out  the  indicated  summation 
using  equation  4.2  gives  the  following  expression  for  the  mgf: 
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Q^(s)  =  exp  [W(exp  (-s)  -1)]  (4.6) 

The  usefulness  of  the  mgf  is  apparent  when  the  term  exp(-ns)  is  expanded 
in  a  Maclaurin  series.  The  mth  term  of  the  expansion  evaluated  at  s  =  0  is 
the  mth  ordinary  moment  of  n.  Using  equation  4.2  and  the  mgf  the  rela- 
tions between  the  moments  of  W  and  n  are  easily  derived.  Several  of  these 
are  listed  here  for  future  reference: 

ordinary  moments 

<  n>  =  W  a) 

<n2>  =  W2+W  b)  (4.7) 

<nin2)  =  W1W2  c) 

central  moments 

<5n)  =  <(n-<n)))  =  0  a) 

<(5n)2)  =  W  b)  (4.8) 

(5ni5n2)  =  0  c) 

The  statistics  of  a  non-homogeneous  P.PP  are  completely  determined  by 
the  rate  density  p(t).  However,  if  P(t)  is  itself  a  stochastic  function  then  the 
statistics  of  the  point  process  depend  on  the  statistics  of  p.  If  p(t)  is  a 
stationary  stochastic  function  then  the  point  process  described  by  it  is  called 
a  doubly  stochastic  Poisson  point  process  (DSP.PP).  The  moments  and 
moment  generating  functions  of  the  DSP.PP  are  found  by  averaging  the 
corresponding  moments  and  mgf  s  of  the  P.PP  over  the  different  realiza- 
tions of  P(t)  or  the  integrated  rate  W.  Thus,  the  moments  of  the  DSP.PP  are 
obtained  by  averaging  equation  4.7  over  W.  They  are  listed  here: 

ordinary  moments  of  DSP.PP 

<n)  =  (W)  a) 

(4.9) 
<n2>  =  <W2)+(W)  b) 
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and,  conversely: 

(W)  =  <n>     ■  a) 

(4.10) 
(W2>  =  <n2>-<n>  b) 

Also: 

<nin2>  =  (WiW2>  (4.11) 

-      variance  of  DSP.PP 

<(5n)2)  =  ((6W)2)+<W)  (4.12) 

These  relations  and  the  results  of  section  III  are  used  to  derive  the  photo- 
electron  counting  statistics  for  the  binary  source  method. 

B.    THE  RATE  OF  PHOTOEMISSIONS 

A  full  quantum  analysis  of  the  photoemission  process  reveals  the 
following  facts  [Ref  25]: 

1.  The  probability  of  transition  from  a  bound  to  an  unbound  state  for 
an  electron  in  a  photosensitive  surface  during  a  time  At  is  proportional  to  At 

and  the  instantaneous  intensity. 

2.  If  the  intensity  is  a  stochastic  function  of  time  then  the  photoemis- 
sion process  is  a  DSP.PP. 

These  conditions  hold  quite  generally  if  the  following  two  auxiliary 
conditions  hold: 

1.  The  time  At,  though  much  greater  than  the  period  of  oscillation  of 

the  EM  wave,  is  much  shorter  than  other  characteristic  time  scales  of  the 
experiment. 

2.  The  bandwidth  of  unbound  (free)  electron  states  is  much  greater 
than  the  EM  bandwidth. 
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Under  these  conditions  then  the  rate  density  of  photoelectric  emissions 
can  be  expressed  as: 

P(t,r)  =  JJ    .     dv  d\  a  (V)  I(r,t,v)  (4.13) 

with  V  the  EM  frequency,  Ad  the  detector  element  area  and  a(v)  the  quan- 
tum efficiency  at  frequency  v. 

This  expression  for  the  rate  density  is  simplified  by  making  the  quasi- 
monochromatic  approximation.  The  quasimonochromatic  approximation 
assumes  that  the  center  EM  frequency  is  much  greater  than  the  EM  band- 
width. For  optical  frequencies  this  condition  is  easily  satisfied  by  using  an 
appropriate  set  of  colored  filters.  However,  the  same  thing  is  accomplished 
by  a  narrowband  detector  response.  If  a  more  accurate  assessment  of  the 
frequency  dependence  is  required  the  stellar  spectrum  is  modeled  by  a 
Planck  blackbody  spectrum  and  the  explicit  frequency  dependence  of  the 
quantum  efficiency  is  used  to  evaluate  the  frequency  integral  of  equation 
4.13.  However,  since  an  analytic  form  for  the  quantum  efficiency  is  not 
readily  available  this  calculation  proceeds  by  invoking  the  quasimonochro- 
matic approximation  and  approximating  the  frequency  integral  as  follows: 

p(r,t)  =  Ava(vjL    d\  I(r,t,v  )  (4.14) 

o     Aj^  o 

with  Vo  the  center  frequency. 

The  spatial  integral  can  similarly  be  approximated  by  assuming  that 
the  area  over  which  the  spatial  fluctuations  are  correlated  is  much  larger 
than  the  detector  element  area.  The  rate  density  then  becomes: 

p(r.,t)  =  A  V  a  (v^)  A^  I(r.,t,v^)  (4.15) 

with  ri  the  geometric  center  of  the  detector  element. 
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The  rate  of  photoemissions  is  given  by  equation  4.3: 

t+T  t+T 

W(r.,t)  =  J  p(t,r.)dt  =  Y  J  I(r.,t)  dt  (4.16) 

t  t 

with  Y  equal  to  Ava(Vo)AD  and  letting  the  frequency  dependence  be  under- 
stood. But,  as  noted  in  section  III,  the  intensity  is  the  sum  of  two  indepen- 
dent sources: 

Kl^.t)  =  I^(r.,t)  +  l2(r.,t)  (4.17) 

Following  Fried  [Ref.  23],  the  intensity  is  expressed  as  a  function  of  the  log 
amplitude  as: 

I.(r.,t)  =  f  (t)  exp  (2x.(r.,t))  (4.18) 

with  Ii°(t)  the  source  intensity;  and,  the  exponential  term  is  the  atmo- 
spheric modulation.  Again,  the  source  term  is  assumed  spatially,  though 
not  temporally,  coherent  over  the  detector  element  area.  Under  these 
assumptions  the  rate  is: 

t+T 

W(r.,t)  =  Y  J  [l°(t)  exp  (2x^(r.,t))  +  I2(t)exp(2x2(r.,t))]dt         (4.19) 
t 

The  stellar  sources  are  assumed  thermal  and  as  is  well  known  thermal 
sources  have  a  coherence  time  that  is  the  inverse  of  the  optical  bandwidth. 
However,  the  correlation  time  associated  with  the  atmospheric  modulation 
is  on  the  order  of  one  millisecond.  This  is  consistent  with  the  assumption  of 
frozen  turbulence  for  times  less  than  or  equal  to  one  millisecond  [Ref.  20]. 
For  any  reasonable  optical  bandwidth  the  source  coherence  time  is  safely 
assumed  to  be  many  orders  of  magnitude  less  than  the  correlation  time  of 
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the  atmospheric  modulation.  Selecting  an  integration  time  that  is  much 
less  than  the  atmospheric  correlation  time  and  much  greater  than  the 
source  coherence  time  the  rate  becomes: 

W(r.,t.)  =  W>.)  exp  (2x.(r.,t.))  +  W°(t.)  exp(2x_(r..t.))  (4.20) 

with  tj  =  t+T/2.  Note  that  the  atmospheric  modulation  is  considered  frozen 
over  the  integration  time. 

Equation  4.20  is  used  to  form  the  second  order  moments  of  the  rate  and 
these  moments  are  related  to  the  second  order  moments  of  the  photo- 
electron  counts. 

C.    COVARIANCE  OF  PHOTOELECTRON  COUNTS 

Forming  the  normalized  spatial  covariance  of  the  integrated  rate  yields: 

<5W(r^.t.)  5W(r_,t.))         <W(r^,t.)  W(r_,t.)) 

LJ ^_J_  =  LJ ^_J —  _  1  (4  21) 

<W(r^,t.))(w(r2,t.))        <W(r^,t.))(W(r2.t.)) 

with  the  angled  brackets  denoting  an  ensemble  average.  Examining  the 
correlation  term  of  equation  4.21  and  noting  that  the  atmospheric  and 
source  modulation  are  independent  then  this  term  becomes: 

(W(r^,t.)  W(r2,t.))  =  (Wj(t.))(exp[2x^(r^,t.)  +  2x^(r2,t.)]) 
+  (W^(t.))<exp[2x2(r^,t.)  +  2x2(r2,t.)])  +  (w°(t.)  W°(t.)) 
•  {(exp[2x^(r^,t.)  +  2x2(r2,t.)])+(exp[2x^(r2,t.)+2x2(r^,t.)])}      (4.22) 

by  using  equation  4.20.  This  relation  is  further  simplified  by  using  the 
relation  between  the  intensity  and  the  log  amplitude.     Following  Fried 
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[Ref.  23],  the  exponential  terms  of  equation  4.22  are  related  to  the  log  ampli- 
tude covariance  by: 

<exp  [2x.(r^,tj^)  +  2x  (r^,tj^)])  =  exp  [4B^    (p)]  (4.23) 

ij 

Since  the  log  amplitude  fluctuations  are  assumed  small  then  this  expres- 
sion is  expanded  in  a  Maclaurin  series  to  yield: 

exp  [4B^   (p)]  =  4B^   (p)  -i-  1  (4.24) 

The  various  log  amplitude  covariances  implied  by  equation  4.24  are  again 
simplified  by  appealing  to  the  symmetry  of  Figure  3  and  equations  3.31  to 
3.34.  The  result  is: 

<W(r^,t.)  W(r2,t.))  =  [(W°(t.)^)  +  (W°(t.)^)]  (4B^(p)  +l) 

+  4<W°(t.)  W°(t.))[B  (p-eh)  +  B  (p+0h)]  (4.25) 

J.  1  ^  1  X  A 

or,  rearranging  terms: 

<W(r^,t.)W(r2,t.))  =  <[W°(t^)+W°(t.)]^)  +  [<W°(t.)^)+<W°(t.)^)]4B^(p) 

+  4(w°(t.  )W2(t. ))  [4B^(p-eh)  +  4B^(p+eh)]  (4.26) 

Noting  that  (exp[2x(r,t)])  equals  unity  by  conservation  of  energy  then  the 
normalizing  term  of  equation  4.21  is  expressed  as: 

<W(r^.t.))(w(r2,t.))  =  [<W°(t.))  +  (W°(t.))]^  (4.27) 

by  using  equation  4.20.  Therefore,  substituting  from  equations  4.26  and  4.27 
in  equation  4.21  yields  the  following  expression  for  the  spatial  covariance  of 
the  rate: 
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<5W(r^,t.)  5W(r2,t.))      <[wj(t.)  +  W^(t.)]^) 
<W(r^,t.))<W(r2,t.))       [<Wj(t.))+<W2°(t.))]' 


[<W°(t.)')+(w^(ty)]  4«(t.)  W°(t.)) 
V4B  (p)+  1^21 


[<W°(t.)Mw°(t.))]        ^        [<wJ(t.)Mw°(t.))]' 

•  [B^(p+0h)  +  B^(p-eh)]  -  1  (4.28) 

As  stated  previously  the  photoelectron  emission  process  is  a  DSP.PP  for 
a  stochastic  intensity.  The  relations  between  the  second  order  moments  of 
W,  the  integrated  rate  and  n,  the  number  of  counts,  are  given  by  equations 
4.9  to  4.12.  Substituting  from  these  relations  in  the  rate  spatial  covariance 
of  equation  4.28  yields: 

<5n(r^,t.)  5n(r^,t.))  ^    ([n^(t.)  -k  n^(t.)]^)        [(n^(t.))  +  (n^(t.))] 
<n(r^,t.))(n(r2,t.))      [(n^(t.))  +  (n^Ct.))]^      [<n^(t.))  +  (n^Ct.))]" 

[(n^itf)  +  (n^itf)  -  <n^(t.))  -  (n^Ct.))] 

+ 2 ^^  ^P^ 

[<n^(t.)')  +  <n2(t.))] 

4<n  (t.))  <n  (t.)) 

^   '         ^   '         [B  (p-0h)  +  B  (p+eh)]  (4.29) 


[(n^(t.))  +  (n2(t.)>] 


2    •-    x^ 


This  is  the  normalized  spatial  covariance  of  photoelectron  counts  for  a 
binary  stellar  source. 

To  proceed  further  note  that  a  single  measurement  or  counting  interval 
T  of  approximately  one  millisecond  probes  a  single  configuration  of  the 
atmospheric  modulation.  This  single  measurement  is  clearly  insufficient 
to   determine   the  modulation   statistics,   and   since   there  is   only  one 
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atmosphere  to  sample  then  the  ensemble  averages  of  equation  4.29  must  be 
replaced  with  a  time  average.  Assuming  ergodicity  holds,  Saleh  [Ref.  28] 
points  out  that  long  time  averages  wash  out  the  stochastic  nature  of  the 
intensity  fluctuations.  To  demonstrate  this  assume  that  the  unmodulated 
source  intensity  is  characteristic  of  linearly  polarized  thermal  light.  Saleh 
derives  the  probability  density  of  photoelectron  counts  in  terms  of  the 
parameter  N: 

P(n)  =  ("^^-^)     (l.M)-'^(l.il)-" 
n  ^^  <n> 

with  N  the  number  of  coherence  times  of  the  source  intensity.  Apparently, 
a  good  determination  of  the  atmospheric  modulation  is  made  only  after 
many  coherence  times  of  the  source  intensity  have  elapsed.  In  the  limit 
that  N  goes  to  infinity  equation  4.30  becomes  the  Poisson  distribution: 

P(n)  =  -^exp(-<n»  (4.31) 

n! 

And  as  is  well  known  [Ref.  26],  the  following  relations  hold  for  the  Poisson 
distribution: 

((6n)^)  =  <n)  (4.32) 

(n^)  =  (n)%(n)  (4.33) 

Therefore  in  the  limit  of  long  (compared  with  one  millisecond)  time  aver- 
ages the  normalized  spatial  covariance  of  counts  is: 
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(6n(r^,t.)5n(r2,t.))      {n)\  <n/  4  (n  )  <n) 
T4B^(p)  + 


•  [B^(p-eh)  +  B^(p+eh)]  (4.34) 

This  expression  indicates  that  the  spatial  covariance  of  counts  is  directly 
dependent  on  the  second  order  statistics  of  the  atmospheric  modulation. 
The  source  modulation  is  averaged  away  due  to  the  extremely  short 
(Av  >10SHz)  coherence  time  of  the  source  radiation. 

Since  the  actual  experiment  depends  on  replacing  ensemble  averages 
with  time  averages  the  statistics  of  the  experimentally  estimated  spatial 
covariance  must  be  explored.  The  next  section  examines  the  counting 
experiment  for  a  single  layer  Ah. 

D.    STATISTICAL  ACCURACY 
The  spatial  covariance  of  counts: 

<5n(r^,t.)  6n(r2,t.))  <n(r_,t.)  n(r.  ,t.)) 

=  ^ 1  (4.35) 

<n(r^,t.))  <n(r2,t.))       <n(r^,t.))  (n(r2,t.)) 

can  be  estimated  by  the  expression: 


A 

G 


with; 


S  =  -. -. -1  (4.36) 

n(r^)  n(r2) 


^  =  fS"""!'*™)  "'■-2-tm>  •>  "(--j)  =  :S-|i  "<^j'tm)  ^.37) 
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N  is  the  number  of  counting  intervals  of  length  T,  and  T  is  expected  to  be 
much  longer  than  the  source  coherence  time  and  less  than  the  atmospheric 
modulation  correlation  time.  The  estimator  g  is  characterized  by  its  mean 
and  variance.  The  mean  or  expected  value  of  g  is  apparently  given  by  equa- 
tion 4.34. 

(g/  = 4B  (p)  +  

•  [B^(p-eh)  +  B^(p+eh)]  (4.38) 

However,  the  variance  of  g  is  a  difficult  quantity  to  derive.  The  variance  of  g 
is  defined  by  the  relation: 

var(g)  =  ((g-gf)  (4.39) 


which  simplifies  to: 


:h-{[-A-:-{-A-)f) 


var  (g)  =  \  It—;—-  -  \ v^x— / J  ;  (4.40) 

n(r^)n(r2)  n(r^)n(r2) 

This  quantity  is  the  same  as  the  variance  of  the  estimated  spatial  correla- 
tion  of  counts.  Saleh  evaluates  this  variance  by  expanding  G  and  hj  about 
their  means  and  retaining  only  the  deviations  from  these  means.  With  the 
additional  assumption  that  N  is  very  large  the  variance  of  g  is  given  by 
equation  7.181  of  Reference  25.  In  the  notation  of  this  calculation,  the  vari- 
ance of  g  is: 
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var(g)  =  var(G)/«n(r^))(n(r2)))  +  A 

A  =  4<g)  +  i:    [ ^ +     ^^^    2   var(n(r  ))J 

J=1'2    (n(r^))<n(r2))<n(r  )>      <n(r  )) 

2  cov  «n(r  j),  (n(r  j)) 

+ (4.41) 

<n(r^))  (nCr^)) 

with  the  indicated  variances  the  full  spatiotemporal  variances  of  the  indi- 
cated quantities.  Unfortunately,  a  complete  evaluation  of  this  expression 
reqiiires  the  time  statistics  of  the  atmospheric  modulation.  The  turbulence 
in  the  atmosphere  is  typically  intermittent  in  intensity  and  the  statistics  of 
intermittency  are  inaccessible  with  current  theoretical  models.  However, 
the  single  scattering  assumption  allows  a  lower  bound  to  be  placed  on  the 
variance  of  g. 

The  time  statistics  of  the  evolving  refractive  index  inhomogeneities  will 
depend  on  moments  of  the  velocity  field  of  higher  order  than  is  represented 
by  the  structure  function  (see  section  II).  This  implies  that  the  temporal 
variance  of  the  spatioangular  correlation  of  log  amplitude  fluctuations  will 
be  a  quantity  that  is  fourth  order  in  the  log  amplitude.  By  the  single 
scattering  approximation  these  moments  are  approximately  zero.  There- 
fore neglecting  the  terms  in  the  variance  of  g  that  are  of  order  g2  and 
higher,  then  equation  4.41  becomes: 

var(g)  >  4^)  + ^ =  4(g)  [l  + '^ -]        (4.42) 

{n(r^))<n(r2))  [(n^Mn^)] 

with  ni,  n2  the  counts  due  to  source  1  and  2  respectively.  The  relative  error 
is  defined  as: 
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1/9 

e  =  100[var(g)]     /<g)  (4.43) 

or, 

am 

A  graph  of  the  relative  error  for  (ni>  =  (n2)  =  n  is  given  in  Figure  8  as  a 
fiinction  of  the  spatioangular  covariance  of  counts  <g>. 

From  Figure  8,  the  relative  error  diverges  rapidly  as  <g>  goes  to  zero. 
Since  observed  values  of  Cn^  are  very  small  (10"12  ^q  lO'l^cml^^)  [Ref.  12] 
then  the  relative  error  induced  due  to  the  detection  process  will  dominate 
the  counting  statistics.  Clearly  some  error  reduction  scheme  is  required  to 
extract  the  spatial  covariance  of  counts. 
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Figure  8.      Relative  Error  Versus  Spatioangular  Covariance 
of  Counts  for  Various  Values  of  Incident  Intensity 
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V.  LEAST  SQUARES  ANALYSIS 

The  analysis  of  the  preceding  section  indicated  the  dominance  of  the 
detection  uncertainty  in  the  determination  of  the  spatioangular  covariance. 
This  large  error  is  ameliorated  by  using  the  redundant  information 
provided  by  the  multiple  detector  elements  of  the  array  (Figure  3).  Least 
squares  analysis  is  employed  to  use  this  redundant  information  to  reduce 
noise. 

Most  of  the  analysis  to  this  point  has  applied  to  a  single  layer  at  a  mean 
altitude  h.  However,  equation  4.38  is  easily  generalized  to  any  layer  with  a 
mean  altitude  of  hi: 

(n  )  +  (n  ) 

4<n^)<n2) 
o    [B  (p.-eh.)  +  B  (p.+0h.)]  (5.1) 

Note  that  the  log  amplitude  covariances  are  now  calculated  at  the  centers  of 
the  detector  elements  pj.    With  the  previous  assumption  of  a  small  detector 

element  size  this  is  a  reasonable  approximation.  Rewriting  the  spatio- 
angular weighting  functions  to  reflect  the  explicit  dependence  on  the 
refractive  index  structure  parameter  yields: 

B^(p.,h.)  =  <C^(h.))  B^(p.,h.)  (5.2) 

and  equation  4.38  is  generalized  to: 
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C^(p.,hp  =  (Clih^))  i      \        '       4B>.,h.)  + 
4(^l)(^2)     r  n\ 

— ^[b  (p.-eh.)  +  B  (p.+eh.)]  f 

=  <C^(h.))R(p.,h.)  (5.3) 

^    n     1  '  J     1 

Using  this  equation,  and  assuming  the  fluctuations  induced  by  non-over- 
lapping layers  are  uncorrelated,  then  spatioangular  covariance  of  counts 
observed  in  the  aperture  due  to  multiple  layers  is  given  by  the  summation: 

CR(p  )  =  E  C^(p  h.)  =  E  <C^(h.))  R(p.,h.)  (5.4) 

•'         h.  "'  h.  "' 

1  1 

This  summation  is  equivalent  to  the  matrix  equation: 

[CRj  =  [R^.]  [(Cl).]  (5.5) 

with  CRk  the  experimentally  measured  spatial  covariance  of  counts,  Cn^j 
the  undetermined  refractive  index  structure  parameter  profile  and  Rjk  the 
weighting  functions  described  above.  The  undetermined  Cn^  are  found  by 
minimizing  the  quantity: 

I  [C\]  -  [Rj^.]  [(C^).]  I    =  minimum  (5.6) 

as  a  function  of  Cn^.  The  solution  to  this  is  well  known  and  is  given  by: 

[cy  =  (Kj]*[V)''Kj]*[CRj  (5.7) 

with  the  superscript  *  denoting  the  Hermitian  conjugate.  This  process  is 
formally  known  as  least  squares  analysis  and  the  Cn^  vector  obtained  is  the 
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least  squares  approximation  of  the  structure  parameter  profile.  A  more 
explicit  expression  for  the  profile  is  obtained  by  considering  the  experimen- 
tally derived  spatioangular  covariance. 

The  experimental  covariance  is  given  by  the  matrix: 


[CR.  J  =  [n(r.)n(r^)  -  l] 


(5.8) 


Assuming  the  angular  separation  of  the  binary  source  is  small  then  the 
expected  values  of  these  matrix  elements  are  a  function  only  of  the  separa- 
tion between  detector  elements.  Therefore,  elements  along  diagonals 
parallel  to  the  main  diagonal  are  expected  to  be  equal  in  the  absence  of 
noise.  Using  this  fact,  the  spatioangular  covariance  of  counts  as  a  function 
of  detector  element  separation  is  formed  by  summing  along  these  diagonals 
as  indicated  here. 


[CRjk]  -^  [CRkl  by 


n(rj^)n(r3^)       nCr^Mrg)  .  .  .  n(rj^)n(r^Q) 
n(r2)n(ri)^  nCrgr  ^^^^^^ 
ii(r3)n(rj)       n(r3)n(r2r-. 


n(r^Q)n(r^) 


^n(r,n)n(r,,) 


10^ 


10^ 


(5.9) 


Since  the  number  of  terms  along  these  diagonals  are  not  equal,  then  the 
elements  of  this  sum  must  be  normalized  by  the  number  of  terms  to 
preserve  the  relationship  of  equation  5.4.  Therefore,  the  elements  of  the 
measured  covariance  vector  are  written: 
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CR,  =(2(10-ln-ml))     Z      [CR       +  CR     ]  (5.10) 

^  n=l,10-k  "™  "^" 

m=n+k 

and,  the  structure  parameter  profile  is  determined  by  solving  the  matrix 
equation 

[|^(\iM   [C>j']  =  [?  ('^kj  CRu)]  (5-11) 

Which  results  from  applying  least  squares  analysis  with  the  parameters 
Cn^(hj)  the  undetermined  coefficients.  The  profile  is  now  determined  since 
the  detector  element  separation  and  altitude  are  connected  by  the  relation 
3.1. 
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VI.    CONCLUSIONS 

This  thesis  has  explored  the  theoretical  basis  for  profiling  the  refi'active 
index  structure  parameter  using  the  spatial  covariance  of  binary  source 
intensity  scintillations.  This  analysis  derived  the  power  spectrum  and 
associated  spatial  covariance  of  intensity  scintillations  caused  by  atmo- 
spheric refractive  turbulence.  The  expected  spatial  covariance  of 
photoelectron  (photon)  counts  and  counting  statistics  were  derived  using 
the  intensity  covariance  functions. 

No  rigorous  analysis  of  the  binary  source  technique  as  implemented  by 
several  French  teams  [Refs.  12-15]  is  currently  available  in  the  open  litera- 
ture; and,  the  analysis  of  this  thesis  reveals  several  interesting  features  not 
previously  noted. 

The  power  spectrum  and  spatial  covariance  used  by  the  French  teams 
are  essentially  correct;  however,  the  error  analysis  of  the  binary  source 
technique  is  incomplete.  A  consideration  of  the  photoelectron  counting 
statistics  indicates  that  very  large  relative  errors  of  200  to  2,000  percent  are 
expected  in  any  determination  of  the  spatial  covariance  function.  This 
large  error  is  not  significantly  improved  by  increasing  the  counting  rate 
above  unity  as  indicated  by  Figure  8.  Previous  experimental  work  [Refs. 
12-15]  was  limited  to  bright  (magnitude  2)  stellar  binaries.  However,  an 
elementary  calculation  of  the  available  photons  per  counting  time 
(approximately  one  msec)  coupled  with  the  relative  insensitivity  of  the 
binary  source  technique  to  increased  counting  rates  indicates  that  stellar 
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binaries  of  magnitudes  three  to  four  are  also  suitable  for  detector  apertures 
of  30  to  60  centimeters. 

The  large  relative  error  expected  in  an  experimental  determination  of 
the  spatial  covariance  necessitates  the  use  of  a  smoothing  algorithm  and 
the  least  squares  technique  was  selected.  The  advantage  of  least  squares 
analysis  is  that  it  accomplishes  the  two-fold  purpose  of  data  smoothing  and 
refractive  index  structure  parameter  profiling.  The  basis  selected  for  the 
least  squares  analysis  is  the  nonorthogonal  set  of  theoretically  derived 
spatial  covariance  functions  for  a  series  of  atmospheric  layers  of  increasing 
attitude.  The  weighting  coefficients  derived  using  the  least  squares  algo- 
rithm are  found  to  be  the  required  refractive  index  structure  parameters  for 
the  atmospheric  layers. 

The  expected  improvement  in  the  relative  error  remains  an  open-ended 
question.  The  error  analysis  depends  on  the  intermittency  of  atmospheric 
turbulence  and  as  indicated  in  Section  IV,  the  error  analysis  in  this  thesis 
constitutes  a  lower  bound  for  the  expected  error.  Clearly,  further  work  on 
the  intermittent  aspects  of  atmospheric  turbulence  is  indicated. 

Despite  the  large  errors  expected  in  the  determination  of  the  spatial 
covariance  function  the  binary  source  technique  has  several  clear 
advantages.  As  mentioned  in  the  Introduction,  the  altitude  resolution  of 
the  binary  source  technique  is  superior  to  that  of  existing  methods.  The 
implementation  of  the  least  squares  algorithm  is  straightforward  and 
results  in  a  relatively  unambiguous  determination  of  the  structure 
parameter  profile.  Modest  aperture  sizes  (30  to  60  cm)  are  adequate  for  the 
detector  system.   The  greatest  disadvantage  of  the  technique  is  the  paucity 
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of  suitable  stellar  binaries.     However,  by  using  stellar  sources  of  up  to 
magnitude  four  this  drawback  becomes  minor. 

Based  on  this  analysis,  experimental  implementation  of  this  technique 
should  proceed,  and  the  experimental  profiles  obtained  should  be  compared 
with  insitu  measurements  or  integrated  profiles  such  as  that  produced  by 
stellar  scintillometers  [Ref  16]. 
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